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Abstract 

Acoustic  sensor  arrays  will  likely  play  a  prominent  role  in  the  U.  S.  Army’s  Objective  Force  for 
applications  such  as  target  detection,  identification,  and  location.  Generally,  these  array  systems  perform 
direction  finding  by  determining  the  wavefront  angle  of  arrival  from  the  phase  differences  across  the 
array.  However,  random  fluctuations  in  the  atmosphere  may  strongly  distort  the  wavefront  of  a  sound 
wave.  The  resulting  random  variations  in  the  wavefront’s  orientation  and  intensity  are  perceived  as 
fluctuations  in  the  apparent  bearing  angles  and  strength  of  the  source.  As  such,  the  error  in  estimating 
the  wavefront’s  angles  of  arrival  (AOAs)  will  increase  as  the  propagation  distance  increases  and/or  the 
strength  of  the  turbulence  increases.  The  Cramer-Rao  lower  bound  (CRLB)  is  the  statistical  lower  bound 
of  the  mean-square  error  between  an  estimator  and  its  actual  value.  The  CRLB  is  calculated  from  the 
probability  likelihood  function  of  the  received  signal.  Previous  formulations  have  modeled  the  probability 
likelihood  function  of  a  sound  wave  propagating  through  atmospheric  turbulence  as  a  zero-mean  complex 
Gaussian  distribution.  The  zero-mean  assumption  implicitly  treats  the  case  of  strong  scattering.  In  this 
paper  the  received  signal  is  modeled  as  a  complex  Gaussian  random  variable  with  a  deterministic  mean, 
thereby  generalizing  the  treatment  to  both  strong  and  weak  scattering.  The  CRLB  of  the  azimuthal 
and  elevational  AOAs  are  calculated  for  both  a  plane-wave  and  a  spherical-wave  propagating  through 
atmospheric  turbulence  with  fluctuations  described  by  a  von  Karman  spectrum.  A  single  monochromatic 
source  is  considered  and  a  line-of-sight  propagation  path  is  assumed. 


Introduction 

Acoustic  sensor  arrays  have  long  been  used  in  underwater  applications  such  as  target  detection,  identifica¬ 
tion,  and  location.  Recently,  there  has  been  much  interest  within  the  U.  S.  Army  to  use  atmospheric  acoustic 
arrays  to  perform  the  same  functions  on  land.  Acoustical  direction-finding  and  tracking  systems  will  likely 
play  a  prominent  role  on  the  future  battlefield,  where  situational  awareness  will  be  a  key  factor  affecting 
the  survivability  of  light-  and  medium-weight  forces.  Acoustic  arrays  are  beneficial  as:  they  are  relatively 
inexpensive  and  small  in  size;  they  operate  passively;  obscurants,  such  as  smoke,  do  not  effect  detection; 
and  a  line-of-sight  propagation  path  is  not  necessary.  Generally,  these  systems  perform  direction  finding  by 
determining  the  wavefront  angle  of  arrival  (AOA)  from  the  phase  differences  across  the  array.  Sound  waves 
are  strongly  affected  by  the  environment,  whether  it  be  oceanic  or  atmospheric.  Both  media  are  random  in 
nature  and  may  strongly  distort  the  wavefront.  The  resulting  random  variations  in  the  wavefront’s  orienta¬ 
tion  and  intensity  are  perceived  as  fluctuations  in  the  apparent  bearing  angles  and  strength  of  the  source. 
As  such,  the  error  in  estimating  the  wavefront’s  AOA  will  increase  as  the  propagation  distance  increases 
and/or  the  intensity  of  the  turbulence  increases.  (These  acoustic  phenomena  are  analogous  to  scintillation 
and  quivering  of  optical  images,  as  are  often  observed  above  a  roadway  on  a  sunny  afternoon.)  The  net 
effects  of  these  distortions  can  have  a  substantial  impact  on  direction- finding  in  both  the  atmosphere  [1,  2] 
and  the  ocean  [3,  4]. 

The  performance  of  a  sensor  array  may  be  quantified  by  calculating  the  mean  square  error  (MSE)  between 
the  estimated  parameter  (such  as  the  AOA)  and  its  actual  value.  The  lower  bound  of  the  MSE  is  the  Cramer- 


Report  Documentation  Page 


Report  Date  Report  Type 

00  Oct  2001  N/A 

Dates  Covered  (from..,  to) 

Title  and  Subtitle 

Contract  Number 

Cramer-Rao  Lower  Bounds  ot  Angle-ol- Arrival  Estimates  for 
Acoustic  Sensor  Arrays  Operating  in  Atmospheric  Turbulence 

Grant  Number 

Program  Element  Number 

Author(s) 

Project  Number 

Task  Number 

Work  Unit  Number 

Performing  Organization  Name(s)  and  Address(es) 

U.S.  Army  Research  Laboratory  ATTN:  AMSRL-CI-EP  2800 
Powder  Mill  Road  Adelphi,  MD  20783-1197 

Performing  Organization  Report  Number 

Sponsoring/Monitoring  Agency  Name(s)  and  Address(es) 

Department  of  the  Army,  CECOM  RDEC  Night  Vision  & 

Sponsor/Monitor’s  Acronym(s) 

Electronic  Sensors  Directorate  AMSEL-RD-NV-D  10221 

Burbeck  Road  Ft.  Belvoir,  VA  22060-5806 

Sponsor/Monitor’s  Report  Number(s) 

Distribution/ Availability  Statement 

Approved  for  public  release,  distribution  unlimited 

Supplementary  Notes 

See  also  ADM201471,  Papers  from  the  Meeting  of  the  MSS  Specialty  Group  on  Battlefield  Acoustic  and  Seismic 
Sensing,  Magnetic  and  Electric  Field  Sensors  (2001)  Held  in  Applied  Physics  Lab,  Johns  Hopkins  Univ,  Laurel,  MD  on 
24-26  Oct  2001.  Volume  2  (Also  includes  1999  and  2000  Meetings),  The  original  document  contains  color  images. 

Abstract 

Subject  Terms 

Report  Classification 

unclassified 

Classification  of  this  page 

unclassified 

Classification  of  Abstract 

unclassified 

Limitation  of  Abstract 

UU 

Number  of  Pages 

20 

Rao  lower  bound  (CRLB),  which  is  calculated  from  the  Fisher  information  (FI)  [5,  6].  There  already  exists 
much  work  in  the  open  literature  that  characterizes  the  CRLB  of  array  processors  in  the  presence  of  noise 
only,  (e.g.,  see  Ottersten  et  al.  [7]  and  the  references  therein).  Recently,  Song  and  Ritcey  [8]  developed 
a  model  that  directly  incorporates  the  effects  of  random  media  on  acoustic  waves  into  the  calculation  of 
the  CRLBs  of  AOA  estimates.  Specifically,  they  considered  propagation  in  an  ocean  channel  with  random 
inhomogeneities  having  a  Gaussian  spatial  correlation.  Using  the  general  framework  of  Song  and  Ritcey, 
Wilson  [9]  calculated  the  performance  bounds  on  AOA  estimates  using  various  correlation  functions  that  are 
representative  of  atmospheric  turbulence. 

Song  and  Ritcey  [8]  and  Wilson  [9]  modeled  the  received  signal  as  a  complex,  zero-mean,  Gaussian 
random  variable.  By  assuming  that  the  signal  of  interest  had  a  zero-mean,  they  implicitly  treated  the  case 
of  waves  strongly  scattered  by  the  turbulence.  A  strong  scattering  event  occurs  when  the  turbulence  is 
sufficiently  strong,  and/or  the  wavefronts  propagate  a  sufficient  distance,  so  that  the  phase  of  the  received 
signal  at  each  sensor  is  completely  randomized.  However,  in  many  practical  problems  the  variance  of  the 
index-of-refraction  fluctuations  is  sufficiently  small,  or  the  propagation  distance  is  sufficiently  short,  that  the 
waves  have  only  been  weakly  scattered  when  they  arrive  at  the  array.  A  weakly  scattered  wave  has  a  mean 
component  significantly  larger  than  the  standard  deviation  of  the  real  and  imaginary  parts  of  the  signal. 

Moreover,  in  the  previous  work  of  Wilson  [9] ,  the  FI  was  calculated  assuming  that  the  only  unknown  was 
the  wavefront  AOA  (azimuth  only).  The  source-receiver  propagation  distance,  turbulence  parameters,  and 
signal-to-noise  ratio  (SNR)  were  implicitly  assumed  to  be  known.  In  a  real  scenario,  this  information  may 
not  be  available.  A  more  relevant  analysis  would  therefore  use  a  multivariate  formalism. 

This  paper  provides  an  analysis  of  the  CRLBs  of  AOA  estimates  for  a  passive  sensor  array  operating  in 
a  turbulent  medium.  The  zero-mean  assumption  is  avoided,  so  that  the  model  is  valid  for  both  strong  and 
weak  scattering.  To  capture  scattering  by  eddies  in  both  the  energy-containing  and  inertial  subranges,  a 
von  Karman  spectrum  is  used  for  the  turbulence  [10].  Furthermore,  the  analysis  from  Ref.  [9]  is  expanded 
to  include  the  propagation  distance,  SNR,  turbulence  parameters,  and  phase  of  the  source  in  the  unknown 
parameter  set.  The  extent  to  which  the  estimates  of  the  AOAs  will  degrade  when  they  are  simultaneously 
estimated  with  other  parameters  is  determined  by  calculating  the  coupling  between  the  estimates.1 

Three-dimensional  (3D)  propagation  is  considered  for  an  incident  plane  wave  and  spherical  wave.  For 
simplicity  a  single  monochromatic  source  and  a  line-of-sight  propagation  path  are  assumed.  A  brief  review 
of  the  phenomology  of  atmospheric  turbulence  is  given  in  Section  1.  The  theoretical  model  is  developed  in 
Section  2  and  the  Cramer- Rao  lower  bound  is  formulated  in  Section  3.  Results  are  given  in  Section  4  and 
conclusions  are  drawn  in  Section  5. 


1  Phenomenology 

The  atmospheric  boundary  layer  generally  refers  to  the  layer  of  air  directly  above  the  earth’s  surface  in 
which  effects,  such  as  heating,  cooling,  and  friction,  are  felt  on  a  short  time  scale  (less  than  a  day).  This  layer 
is  usually  one-half  a  kilometer  to  several  kilometers  thick.  Fluctuations  in  momentum,  heat,  or  matter  may 
result  in  turbulent  motions  that  have  scales  on  the  order  of  the  depth  of  the  boundary  layer  or  less.  Some 
characteristics  of  turbulent  flows  are:  they  have  purely  rotational  motion;  they  are  three-dimensional;  and 
they  are  dissipative  and  thus  require  energy  to  be  sustained.  When  a  tangential  force  exists  at  a  boundary,  a 
shear  turbulence  occurs.  Buoyant  turbulent  motion  may  result,  e.g.,  from  the  thermal  heating  of  the  ground 
which  induces  an  upward  motion  of  the  air.  Turbulent  motion  may  be  characterized  by  an  inner  scale  length 
and  an  outer  scale  length  (C).  In  between  the  inner  and  outer  scales  lies  the  inertial  subrange,  characterized 
by  a  cascade  of  energy  to  progressively  smaller  scales.  (See,  e.g.,  Reference  [11].) 

1  The  theoretical  treatment  in  this  paper  applies  to  the  AOAs,  which  represent  the  orientation  of  the  wavefront  normal  when 
the  sound  reaches  the  array.  For  propagation  in  the  atmosphere,  the  average  horizontal  (azimuthal)  AOA  is  usually  very  close 
to  the  actual  angle  of  bearing  (AOB)  of  the  source,  thereby  making  acoustic  arrays  well  suited  to  determining  the  horizontal 
position  of  a  source.  However,  for  the  vertical  (elevational)  AOA,  the  situation  is  quite  different.  Atmospheric  refraction  may 
bend  the  soundwaves  upward  or  downward,  thus  interfering  with  the  ability  to  determine  the  AOB.  As  such,  most  existing 
acoustic  ground  sensors  do  not  attempt  to  determine  the  elevation  of  a  near-ground  source. 


Figure  1:  Turbulence-induced  distortion  of  acoustic  wavefronts  impinging  on  a  sensor  array. 


As  an  acoustic  wavefront  propagates  from  the  source  to  the  sensor  array,  it  is  distorted  by  random 
fluctuations  in  the  atmosphere.  A  coherence  loss  results  when  the  phase  and  amplitude  relationship  between 
the  sensors  becomes  random.  Let  the  separation  between  the  sensors  be  d ,  and  let  £  be  the  characteristic 
length  of  the  largest  turbulent  eddies.  For  d  <^  £  the  wavefronts  will  be  smooth  and  nearly  planar  when 
arriving  at  the  sensor  array.  But  if  d  >  £,  the  wavefronts  will  have  a  rough  and  random  appearance  when 
observed  on  the  scale  of  the  array.  These  cases  are  depicted  in  Fig.  1.  In  both  cases,  there  is  a  coherence 
loss  between  the  sensor  signals. 

Knowing  £  is  essential  to  designing  a  sensor  array.  Typical  values  of  £  during  the  daytime  range  between 
50  m  (for  strong  winds  and  little  solar  heating  of  ground)  and  500  m  (for  sunny  afternoons  with  light  winds) . 
Most  sensor  arrays  employed  by  the  Army  have  an  aperature  of  about  lm.  Therefore,  the  case  d  <C  £ 
generally  applies.  It  is  also  generally  desirable  to  space  array  elements  less  than  one-half  a  wavelength  apart 
to  prevent  spurious  source  images  (e.g.,  spatial  aliasing  effects). 


2  Theoretical  Model 

Calculation  of  the  FI  requires  a  priori  knowledge  of  the  probability  density  function  (PDF)  of  the  received 
signal.  In  this  section  a  theoretical  model  that  incorporates  the  effects  of  turbulence  on  the  source  signal 
is  developed  to  describe  the  PDF  of  the  received  signal.  Let  us  define  the  notation  that  shall  be  used 
throughout  this  paper:  [•]*  denotes  the  complex  conjugate,  [-]T  the  transpose,  [•]'  the  Hermitian  adjoint 
(complex  conjugate  transpose),  and  (•)  the  ensemble  average. 

Consider  an  acoustic  array  with  N  sensors.  We  assume  that  the  signal  at  each  sensor  results  from:  (1)  the 
wave  that  has  propagated  from  the  source  of  interest  with  <f>  and  9  as  the  azimuthal  and  elevational  angles 
of  arrival,  respectively,  and  (2)  random  noise.  Let  p (cf>,9,t)  and  n(t)  be  the  time-varying  complex  envelops 
of  the  two  contributions,  respectively.  These  column  vectors  have  N  elements,  one  element  corresponding  to 
each  sensor.  The  source  contribution  is  time  dependent  because  of  the  random  turbulent  effects.  The  noise, 
which  is  also  time  dependent,  may  result  from  wind  noise  or  other  competing  sources.  The  total  received 
signal  is 

s((/>,9,t)  =  p((f>,9,t)  +  n(i) .  (1) 

Let  us  assume  that  the  source  and  the  noise  signals  are  uncorrelated  and  that  the  noise  signals  at  the 
sensors  are  mutually  uncorrelated  with  equal  variance.  Let  us  further  assume  that  the  noise  component 
has  a  Gaussian  distribution  with  a  zero  mean  and  variance  cr^.  Although  exact  solutions  for  the  pressure 
field  of  the  source  and  its  PDF  are  unknown,  solutions  to  its  first  and  second  moments  can  be  found  in  the 
literature.  We  therefore  approximate  that  p  has  a  Gaussian  distribution  with  mean  p  and  covariance  matrix 
C„,  whose  elements  are 


Mi  =  (pi)  and  [Cp]..  =  ( PiP* )  -  MiMj  • 


(2) 


We  use  the  results  in  the  open  literature  for  acoustic  wave  propagation  in  a  random  medium  to  determine  /q 
and  (piPj)  as  discussed  in  the  following  section.  The  total  signal  s  is  thus  Gaussian  distributed  with  mean 
p,  and  covariance 

C  =  Cp  +  a2nIN.  (3) 

This  signal  model,  in  which  the  real  and  imaginary  parts  are  modeled  as  Gaussian  random  variables  with 
equal  variances,  is  reasonable  for  strong  scattering,  or  for  weak  scattering,  in  the  presence  of  strong  diffraction 
(the  Rytov  extension  region).  It  is  less  well  suited  to  situations  where  both  scattering  and  diffraction  are 
weak  (geometric  acoustics),  in  which  case  the  phase  variance  dominates  the  signal  behavior  [12]. 


2.1  First  and  Second  Moments:  Normal  Incidence 


The  pressure  field  associated  with  a  sound  wave  propagating  in  a  moving  medium  is  characterized  by  a 
closed  set  of  fluid  dynamic  equations.  The  small  angle  parabolic  equation  method  may  be  used  to  reduce  this 
set  of  equations  to  a  single  wave  equation.  The  Markov  approximation,  which  assumes  that  the  turbulence 
field  has  vanishing  correlation  in  the  direction  of  propagation,  may  then  be  used  to  obtain  the  statistical 
moments  of  the  sound  field  in  closed  form.  These  approximations  are  valid  in  far  field,  for  small  scattering 
angles,  and  for  C  A  >  i,  where  A  is  the  wavelength,  and  C  and  £  are,  respectively,  the  outer  (integral)  and 
inner  (Kolmogorov)  length  scales  of  the  turbulence.  We  use  the  solution  for  the  first  and  second  moments  of 
the  pressure  field  as  given  by  Ostashev  [10],  who  generalized  the  results  of  [13]  and  [14]  to  include  fluctuations 
in  the  medium  velocity.  The  solution  for  the  second  moment  is,  however,  strictly  valid  for  normal  incidence 
across  two  sensors.  These  results  are  outlined  in  the  following  paragraphs. 

By  using  the  small  angle  parabolic  and  Markov  approximations,  one  finds  that  the  first  and  second 
moments  of  the  pressure  field  undergo  an  exponential  attenuation.  Consider  a  sound  wave  that  is  propagating 
along  the  cc-axis  with  wave  number  k.  k  =  Let  r  =  [a;,  y,  z]  be  the  observation  point  and  rx  =  [0,  y,  z] 

be  its  component  transverse  to  k.  For  x  r^,  the  first  moment  at  r  is 

(  p(r)  >=  Pn  (r)  e~lx  ,  (4) 

where  7  is  the  extinction  coefficient  for  the  first  moment  and  pH( r)  is  the  sound  pressure  field  in  the  absence 
of  random  inhomogeneities.  For  a  plane  wave 

PH(r)  =  Po  ei(k  r+x^  ,  (5) 


where  po  is  the  pressure  amplitude  (a  real-valued  constant)  and  %  is  the  phase  of  the  source.  Here  the 
propagation  is  strictly  along  the  x-axis,  therefore  pH( r)  =  po  e^kx+x\  For  a  spherical  wave 

Ph( r)  =  Po  el(fcr+x) ,  where  p0  =  ,  (6) 

r 

and  where  A  is  the  pressure  at  r  =  r0  (A  and  tq  are  real- valued  constants).  The  parabolic  approximation 
for  a  spherical  wave  propagating  in  free  space  (used  in  the  solution  to  the  wave  equation)  is 


exp  (i  kr) 
r 


1 

x 


exp  (\kx)  exp 


(7) 


It  is  the  first  order  approximation  in  /x. 

The  first  moment  represents  the  unscattered  (deterministic)  part  of  the  wavefield.  We  define  the  satura¬ 
tion  parameter  parameter  f l  as 

n  =  1  —  e~2lx  .  (8) 

When  jx  <C  1,  the  saturation  parameter  is  close  to  0,  scattering  is  small,  and  the  signal  at  a  single  sensor 
exhibits  little  variability.  When  ryx  1,  H  ss  1,  and  the  scattered  part  of  the  field  dominates. 

Consider  now  two  points  near  the  x-axis,  iq  =  [x,  yi,  Zi]  and  r-2  =  [x,  2/2,  z 2],  where  x  r^,  ■  The 

second  moment  is 

(  p  (n)  p *  (r2)  )  =  pH  (n)  p*H  (r2)  ,  (9) 


z 


Figure  2:  Coordinate  system.  The  closed  circles  represent  the  sensors.  For  a  plane  wave,  the  open  circle 
represents  the  point  in  the  plane  of  the  source  at  which  r  is  normal  to  the  wavefront.  For  a  spherical  wave, 
the  open  circle  represents  the  source.  The  azimuthal  and  elevational  AOAs,  <fi  and  9,  are  defined  with  respect 

to  r . 


where  p  =  ri  —  r2  is  the  sensor  separation  vector  (transverse  to  the  propagation  direction),  a  is  the  extinction 
coefficient  for  the  second  moment,  and  pH  (17)  is  given  by  Eq.  (5)  for  a  plane  wave  and  Eq.  (6)  for  a  spherical 
wave.  [Again  the  parabolic  approximation  for  a  spherical  wave  propagating  in  free  space,  Eq.  (7),  was  used 
to  obtain  this  result.  It  is,  however,  more  convenient  to  use  the  expression  in  Eq.  (6).  Provided  we  stay 
within  the  constraints  of  the  approximation,  it  does  not  matter  which  expression  we  use.]  The  extinction 
coefficients  for  the  first  and  second  moments  are  related  by 

7  =  a(oo) /2  .  (10) 

2.2  First  and  Second  Moments:  Oblique  Incidence 

We  now  wish  to  derive  expressions  for  the  first  and  second  moments  for  nominally  oblique  incidence. 
Consider  a  sound  wave  that  is  propagating  with  wave  number  k,  where  k  is  no  longer  parallel  to  the  a:-axis. 
Let  [■]_L  and  [■]"  denote  the  components  of  a  vector  transverse  and  parallel  to  k,  respectively.  For  both  a 
plane  wave  and  a  spherical  wave,  r  is  the  propagation  distance  of  the  wavefront  to  the  center  of  the  array. 
However,  the  definition  of  the  vector  r  depends  upon  the  case  considered:  for  a  plane  wave,  r  is  the  vector 
from  the  center  of  the  array  normal  to  the  plane  of  the  source;  for  a  spherical  wave,  r  is  the  vector  from  the 
center  of  the  array  to  the  source.  A  list  of  vectors  and  their  definitions  is  given  in  Table  1.  For  both  a  plane 
wave  and  a  spherical  wave,  the  azimuthal  and  elevational  AOAs,  <f>  and  9 ,  are  measured  with  respect  to  the 
center  of  the  array,  so  that  r  =  [cos</>cos0,  sin^cosf?,  sin0]T.  An  illustration  is  given  in  Fig.  2. 

For  this  analysis  we  take  the  first  moment  at  the  ith  sensor  to  be 


Plane  Wave  : 

Pi  »  poei(k-ri+*}e-^ , 

Po  =  constant 

(11) 

Spherical  Wave  : 

pt=Poe^kri+x)e-^r , 

Ar0 

Po  ~  - • 

(12) 

For  both  the  plane  wave  and  spherical  wave,  the  phase  of  the  first  moment  is  not  modified,  but  the  magnitude 
of  each  is.  In  other  words,  we  assume  that  the  largest  source  of  variation  of  the  first  moment  between  the 
center  of  the  array  and  the  *th  sensor  is  due  to  the  phase  not  the  attenuation.  These  approximations  are 
valid  provided  that  for  every  i  and  j 


Plane  Wave  : 

e~1Ti  «  e_7rJ  «  e“7r 

(13) 

Spherical  Wave  : 

e_7n  ~  e~1Ti  ~  e~1T 

(14) 

It  follows  that  the  saturation  parameter,  for  both  a  plane  wave  and  a  spherical  wave,  is  now 


n»l-  e~2ir . 


(15) 


Plane  Wave 

Spherical  Wave 

r  =  [x,  y,  z] r 

Vector  from  the  center  of  the  array  nor¬ 
mal  to  the  plane  of  the  source 

Vector  from  the  center  of  array  to  the 
source 

x  =  cos  q i  cos  6 

■y  =  sin  (j)  cos  9 

Components  of  r  in  terms  of  the  azimuthal  and  elevational  angles  of  arrival  <f>  and  9 

z  =  sin  9 

r  =  |r| 

Propagation  distance  of  the  wavefront  to  the  center  of  the  array 

r  =  r/r 

Unit  vector  in 

direction  of  r 

G-p 

N/A 

Radial  orthogonal  unit  vector  in  spheri¬ 
cal  coordinates 

k 

Wave  number.  |k|  =  k  =  27t/A,  where  A  is  the  wavelength 

k  =  —  fcr 

k  =  ker 

r'  =  [x't,  yl,  0']T 

Vector  from  the  center  of  the  array  to  the  itli  sensor 

r  j  =  r  +  r' 

As  defined 

Vector  from  the  source  to  the  zth  sensor 

n  =  |r*| 

As  defined 

Propagation  distance  of  the  wavefront  to 
the  zth  sensor 

Sh 

II 

Propagation  distance  of  the  wavefront  to 
the  zth  sensor 

rf  =  ri 

X 

II 

Magnitude  of  the  component  of  r  j  that  is 
transverse  to  direction  of  wave  propaga¬ 
tion 

o 

II 

- l-<s> 

o  =  r'  —  r' 

Hij  A  j 

Vector  between  the  zth  and  jth  sensors 

Pij  =  \Pij\ 

Distance  between  the  zth  and  jth  sensors 

Pij  =  \Pij  x  k| 

Magnitude  of  the  component  of  ptJ  that  is  transverse  to  the  direction  of  wave  prop¬ 
agation 

Table  1:  Definitions  of  vectors  and  their  magnitudes.  Physical  descriptions  are  provided  when  possible. 
Refer  to  Fig.  2. 

Strictly  speaking,  Eq.  (9)  for  the  second  moment  is  valid  for  normal  incidence  across  sensors  i  and  j,  i.e., 
when  r|  =  rj  and  hence  and  hence  pfj  =  py .  For  oblique  incidence  we  approximate 

Plane  Wave  :  {PiP*j)  ~  plelk'l'1'i~rj'>e~a('pij^r  ,  po  =  constant  (16) 

Spherical  Wave  :  (PiPj)  ~  plelk(-ri~rj'>e~a('pij'>r ,  p0  ~  — —  .  (17) 

We  are  thereby  assuming  that  for  every  z  and  j, 

a(Pij)  ~  a(Pij)  (18) 

and 


Plane  Wave  : 


e -<x(pij)r}  ~  e~a{pij)rll  ~  g -a(Pij)r 
e-a(pij)ri  e—Oi{Pij  )rj  e-a(Pij)r 

Spherical  Wave  :  - «  - «  - - - 

nrj  nrj  rz 

A  careful  and  consistent  treatment  of  the  magnitudes  of  the  moments  is  necessary  to  ensure  a  nonsingular 
covariance  matrix.  The  phases  of  the  first  and  second  moments  must  also  be  treated  consistently  to  ensure 
that  the  covariance  matrix  is  nonsingular.  Because  of  these  approximations,  we  restrict  our  investigation  to 
nominally  normal  incidence  to  a  planar  array  for  3D  propagation  [a  linear  array  for  two-dimensional  (2D) 
propagation] . 


(19) 

(20) 


2.3  Covariance  Matrix  and  Mean 

We  write  the  first  and  second  moments  in  the  form 

m  =  Aei$i  and  <  PiP*  )=  p2e~a^re^ 

where  we  define2 


P  =  p0e  ir 


and 


=  •!>,  -  •!>,  • 


For  a  plane  wave  po  is  as  defined  in  Eq.  (5)  and 

<F,;  =  k  •  r?;  +  y  =  k  (r  +  x\  cos  </>  cos  9  +  j/'  sin  (j>  cos  9  +  z\  sin  9)  +  x ; 

and  for  a  spherical  wave  po  is  as  defined  in  Eq.  (12)  and 

r  2  i1/2 

=  h'i  +  x  =  k  r2  +r'i  +2 r  (x'i  cos  (j)  cos  9  +  y\  sin  (f>  cos  9  +  z!i  sin  9) 

Thus  the  elements  of  the  total  covariance  matrix  [Eq.  (3)]  are 


„2  -  2  ,  2 

(-'a  —  Po  p  T  cy 


and 


Ca  = 


p2  e~oc(pij)  r  _ 


9' 


A&i 


(21) 

(22) 

(23) 

(24) 

(25) 


The  mutual  coherence  function  (MCF)  between  the  itli  and  jth  sensors  is  defined  to  be  the  positive 
square  root  of 

( PiPj  )  < P*P j ) 


r2  = 

13  ( PiPi )  (PjPj) 


Thus  for  both  the  plane-wave  and  spherical-wave  treatments  here, 

r«  -  I  {PiP*)  \/Po  =  e-Q(^')r. 

The  minimum  value  of  the  MCF  for  p,j  — >  oo  is  rm;n  =  e^2jr. 


(26) 


(27) 


An  advantage  to  this  formulation  of  the  first  and  second  moments  for  the  plane-wave  treatment  is  that 
the  resulting  covariance  matrix  Cp  and  mean  p,  may  be  written  in  terms  of  the  MCF  matrix  T,  a  steering 
vector  s,  and  a  steering  matrix  S.  The  steering  vector  is  defined  to  be 


s  =  [eik  ri,  eik  r2, 


Dik'rjvl  T 


(28) 


The  steering  matrix  represents  the  phase  delay  between  the  sensors  due  solely  to  the  difference  in  propagation 
length.  As  S  =  S  £g>  where  <S>  is  the  (right)  Kronecker  product, 


Sij  =  exp  (ik  •  ptJ)  . 


(29) 


2Note  that  we  use  the  convention  that  fii  is  the  2th  component  of  the  vector  /x,  (i.e.,  the  value  of  the  first  moment  of  the 
sound  field  at  the  ith  sensor),  and  is  given  by  m  =  jie1®1 .  In  this  manner,  |/Xj|  =  fi  and  |/x|  =  fi  =  y/np,. 


Therefore  for  an  incident  plane  wave 


fj,  =  PuT'n(‘fn  s  and  Cp  =  p20T  ©  S  -  PoFmin<S ,  (30) 

where  ©  is  the  Hadamard  product  (element-by-element  multiplication). 

We  may  also  write  the  moments  for  a  spherical  wave  in  the  same  form  as  Eq.  (30),  but  s  and  S  are 
defined  from  the  phases: 

s  =  [eikri ,  eikr2 ,  •  •  •  ,  eiferjv]T  and  S  =  s  ©  -►  Sxj  =  exp  [ijfc  (r<  -  r,)]  .  (31) 

Though  the  meaning  of  a  steering  vector  and  array  is  lost,  these  mathematical  expressions  may  be  useful  for 
computational  purposes. 


2.4  von  Karman  Turbulence  Model 

The  extinction  coefficients  depend  on  the  structure  of  the  random  medium.  For  a  plane  wave 

a(p)  =  2t rfc2  [/( 0)  -  f{p)\  , 


(32) 


and  for  a  spherical  wave 

a(p)  =  2t rfc2  f  [/( 0)  -  f(pu)]  du  ,  (33) 

Jo 

where  /  is  the  2D  (or  projected)  correlation  function  for  the  sound-speed  fluctuations  [10,  15,  16,  17].  For 
most  random  media,  including  a  turbulent  atmosphere,  a(p)  initially  increases  monotonically  with  increasing 
p ,  but  when  p  exceeds  C,  a(p)  asymptotically  approaches  a  constant  value  [15].  Since  f(p)  — >  0  in  the  limit 
p  — >  oo,  this  constant  value  is  simply  2y,  given  by 

27  =  27rfc2/(0)  =  2 <;2k2C  ,  (34) 


where  o2  is  the  index-of-refraction  variance.  Hence  the  second  moment  will  initially  decrease  with  increasing 
p ,  but  will  eventually  “saturate”  at  a  fixed  minimum  value. 

Typical  acoustic  sensor  arrays  employed  by  the  Army  have  a  sensor  spacing  larger  than  the  height  of  the 
array  from  ground.  As  such,  the  performance  of  these  arrays  is  affected  by  the  large  eddies  of  the  energy- 
containing  (or  source)  subranges  of  the  turbulence  spectrum.  (By  contrast,  the  performance  of  optical 
systems  is  dependent  primarily  upon  the  smaller-scale  eddies  in  the  inertial  and  dissipation  subranges.) 
The  isotropic,  homogeneous  von  Karman  turbulence  model  describes  the  inertial  subrange  of  the  turbulence 
spectrum  more  realistically  than  the  commonly  used  Gaussian  models,  and  it  still  behaves  reasonably  in  the 
energy-containing  subrange.  The  von  Karman  form  for  the  2D  correlation  function  is  dependent  upon  the 
source  of  the  sound  speed  fluctuations:  a  scalar  field  is  induced  by  temperature  or  humidity  fluctuations  and 
a  vector  field  is  induced  by  wind  velocity  fluctuations.  The  2D  correlation  functions  for  a  scalar  field  fs  and 
a  vector  field  fv ,  may  be  written  in  the  form  (see  Eq.  (49)  in  Ref.  [17]  and  Eq.  (7.112)  in  Ref.  [10]) 


fs(p,<;2,l ) 


2c2/ 

(i/s) 

2c2/ 

^r(i/3) 


(35) 

(36) 


where  /  =  r(l/3)£/  [v/7rr(5/6)]  is  a  characteristic  length  scale,  r(cc)  is  the  Gamma  function,  and  Kv(x)  is 
the  modified  Bessel  function  of  order  v.  The  equation  for  /,  in  this  paper  differs  somewhat  from  that  of 
Ref.  [9]  due  to  the  manner  in  which  the  scalar  energy  spectrum  was  defined  in  that  paper.  The  equation  in 
this  paper  is  consistent  with  the  more  standard  definition  found  in  Ref.  [10]. 

The  MCF  for  an  incident  plane  wave  is  plotted  in  Fig.  3  as  a  function  of  the  index-of-refraction  variance 
C2  and  the  characterstic  length  scale  normalized  by  the  wavelength  1/ A  for  both  a  scalar  and  a  vector 


von  Karman  spectrum.  In  presenting  the  results,  it  is  natural  to  use  normalized  length  scales,  (e.g.,  r/A, 
d/X,  etc.),  as  then  the  coherence  has  no  explicit  wavelength  dependence.  In  (a)  and  (b)  the  MCFs  are 
calculated  for  p/X  =  0.5  and  r/A  =  500.  [For  example,  a  tank  has  a  frequency  of  roughly  200  Hz.  Thus, 
as  the  wavelength  is  the  average  sound  speed  (about  340  m/s)  divided  by  the  frequency,  A  is  roughly  2  m. 
Therefore  if  the  tank  is  at  a  range  of  one  kilometer,  r/A  ~  500.]  The  coherence  for  both  spectra  decreases 
significantly  in  the  regions  where  the  index-of-refraction  variance  is  large,  c2  ~  10~4,  and  the  normalized 
length  scale  is  small,  10  <  l/X  <  1.  In  (c)  and  (d)  the  same  is  calculated  but  for  p/X  =  3/\/2.  The  larger 
sensor  separation  leads  to  a  more  rapid  decrease  in  the  MCFs  as  functions  of  the  turbulence  parameters.  For 
both  sensor  separations  the  MCF  for  the  vector  spectrum  is  more  sensitive  to  the  changes  in  the  turbulence 
parameters,  and  its  minimum  with  respect  to  the  turbulence  parameters  (for  a  fixed  finite  sensor  separation 
and  normalized  propagation  distance)  is  smaller  than  that  for  the  scalar  spectrum.  [We  note  that  a(pij)r 
is  dependent  upon  the  product  of  <;2r/X.  Therefore,  one  may  view  the  rr-axis  in  Fig.  3  as  a  change  in  the 
product  of  c2r/A.]  In  Fig.  4  the  same  is  plotted  but  for  an  incident  spherical  wave.  The  function  rm;n 
(the  minimum  value  of  the  MCF  as  a  function  of  sensor  separation  {p  =  oo)  for  fixed  propagation  distance 
and  fixed  turbulence  parameters)  is  the  same  for  both  a  scalar  and  a  vector  spectrum  and  for  both  a  plane 
wave  and  spherical  wave.  It  is  plotted  in  Fig.  5.  Even  though  its  value  is  only  dependent  upon  the  product 
<;2rl/ A2,  it  is  plotted  versus  the  turbulence  parameters  at  r/A  =  500  for  ease  of  comparison  with  Figs.  3  and 
4. 


3  Formulation 


3.1  Cramer-Rao  Lower  Bound 

Consider  the  vector  parameter  ©  =  [0i,  @2,  ...,  ©At]T  that  we  wish  to  estimate.  The  minimum  MSE  of 
an  unbiased  estimator  ©  about  its  actual  value  ©  may  be  calculated  from  the  Cramer-Rao  theorem  [5,  6], 
which  gives 

(  (0,,-ec)2  )>  [J_1(0)]^  ,  (37) 

where  J(@)  is  the  N  x  N  FI  matrix.  (An  estimator  is  said  to  be  unbiased  if  and  only  if  (©)  =  ©.)  The  FI 
is  related  to  the  probability  likelihood  p(x;  ©)  (PDF  of  x  with  ©  as  a  parameter)  by 


[J(©)W  = 


*92lnp(x;  ©)  \ 

90A90„  / ’ 


(38) 


where  the  expectation  value  is  taken  with  respect  to  p(x;  ©)  and  the  derivatives  are  evaluated  at  the  true 
value  of  ©. 


The  likelihood  function  for  real  parameters  of  complex  Gaussian  PDF  with  covariance  matrix  Cx(©) 
and  mean  p(&)  may  be  written  [5] 


p(X;0)  7r^det[Cx(0)] 

Its  corresponding  FI  is  [5] 


exp|-[x-/x(©)]tCx1(©)[x-^(©)]|  . 


[J(©)]a„  =  tr  C 


9C  .  9C 

90a  d®v 


29? 


(dp)_  ±  dp 1 

V  90  a  90  y 


(39) 


(40) 


where  the  functional  dependence  has  been  suppressed.  If  there  are  M  independent  and  identically  distributed 
data  sets,  then  the  likelihood  function  will  be  the  product  of  M  identical  distribution  functions,  and  hence 
the  FI  is  M  times  the  quantity  given  in  Eq.  (40).  Suppose  there  are  N  elements  in  the  sensor  array.  Then 
C  is  an  N  x  N  matrix,  and  p  is  a  column  vector  of  length  N.  Let  us  use  the  convention  that  the  subscripts 
A,  v  £  [1,  . . . ,  J\f]  are  the  indices  on  the  parameters  and  that  i,  j  £  [1,  . . . ,  N]  are  the  indices  on  the  elements 
of  the  sensor  array. 


Let  us  define  oy  =  ^/[J-1]^ .  We  loosely  refer  to  either  a  or  a2  as  the  CRLB,  as  the  meaning  should  be 


(b)  p/X=0.5 


(a)  p/X=0.5 


Figure  3:  Coherence  for  an  incident  plane  wave:  (a)  and  (c)  are  for  a  scalar  von  Karman  spectrum;  (b)  and 
(d)  are  for  a  vector  von  Karman  spectrum.  All  calculations  are  for  r/A  =  500. 


(a)  p/X=0.5 


(b)  p/X=0.5 


Figure  4:  Coherence  for  an  incident  spherical  wave:  (a)  and  (c)  are  for  a  scalar  von  Karman  spectrum;  (b) 
and  (d)  are  for  a  vector  von  Karman  spectrum.  All  calculations  are  for  r/A  =  500. 


evident  from  the  units  involved.  The  minimum  value  of  <r2  is  cr2o  =  1  /  JVV1  i.e.,  the  CRLB  when  0„  is  the 
only  unknown.  As  the  number  of  unknowns  increases,  cr2  will  increase. 


For  example,  suppose  that  there  are  two  unknowns.  The  FI  is 


J  = 


J 11  J 12 
J 12  J22 


(41) 


as  it  is  symmetric  for  real  unknown  parameters.  For  A  and  v  cyclic  (i.e.,  if  A  =  1  then  v  =  2,  etc.) 

=  (Jxx  -  JlJJwY1  =  °l0  (1  -  C12)-1  (42) 


where 


C12  = 


T2 

J12 


=  1  — ^  =  1  — 0<Cl2<l. 


(43) 


J11J22 

Only  if  J12  =  0  does  er2  =  cr2o,  and  the  estimates  of  0i  and  02  are  said  to  be  uncoupled.  As  C12  increases, 
Y  increases  from  its  minimum  value  of  <r2o,  and  a  degradation  of  the  estimates  of  0i  and  02  results.  The 
quantity  £12  thus  provides  a  measure  of  the  strength  of  the  coupling  between,  and  hence  degradation  of, 
the  estimates  of  0i  and  02:  if  £i2  =  0,  the  estimates  are  uncoupled  and  the  CRLBs  retain  their  minimum 
values;  if  £12  <C  1,  the  estimates  of  6\  and  02  are  weakly  coupled  and  the  CRLBs  increase  only  slightly;  and 
if  C12  =  1,  the  estimates  are  fully  coupled,  the  CRLBs  are  infinite,  and  hence  neither  0!  nor  02  can  be 
estimated.  It  is  therefore  advantageous  to  determine  the  conditions  underwhich  the  estimates  of  0i  and  02 
will  decouple.  If  there  are  more  than  two  coupled  parameter  estimates,  we  define  the  coupling  between  the 
Ath  and  iAh  parameter  estimates  to  be 


Cxu  — 


J2 

JXv 


JxxJvv 


(44) 


In  this  way,  we  have  a  measure  of  the  coupling  strength  between  any  two  given  parameters. 


3.2  Fisher  Information  of  Theoretical  Model 

The  dependence  of  the  FI  on  all  the  unknown  parameters  except  the  signal-to-noise  ratio  has  clearly 
been  established.  The  SNR  is  related  to  the  noise  variance  by  SNR  =  pg/cr2.  h  is  often  useful  to  express 
the  SNR  in  decibels  SNR(jn,  SNR  =  ioSNRdB/10.  For  the  plane  wave  case,  po  is  a  constant,  and  we  treat 
the  SNR  as  the  unknown.  In  this  way,  the  exact  value  of  po  is  not  needed  for  the  calculation,  as  the  FI 
may  be  renormalized  by  pg.  For  a  spherical  wave  po  is  dependent  upon  r,  therefore  we  consider  SNR0,  the 
signal-to-noise  ratio  at  a  distance  TZq,  as  the  unknown.  Then  o\  =  (Atro)2/(SNR07?.Q).  By  renormalizing  the 
FI  by  (*4r0)2,  the  explicit  value  of  Ar$  is  not  needed. 

The  FI  may  be  calculated  from  Eq.  (40)  for  those  parameters  we  wish  to  consider  as  unknowns:  </>,  0,  A 
l,  <r2,  or  SNR.  For  a  deterministic  mean,  \  must  be  treated  as  an  unknown  parameter  [22,  23,  24].  While  the 


derivatives  of  the  covariance  matrix  and  mean  are  straight  forward,  the  derivatives  with  respect  to  turbulence 
parameter  l  are  tedious.  For  brevity,  none  of  the  derivatives  are  presented  here.  We  may  interchangeably 
refer  to  the  CRLB  of  the  elevation  and  zenith,  for  as  they  are  related  by  a  linear  transformation,  their  CRLBs 
are  the  same.  (Linear  transformations  preserve  the  efficiency  of  an  estimator  [5].) 


3.3  No  Turbulence 

Let  us  begin  by  examining  the  case  of  no  turbulence.  In  the  absence  of  turbulence 

H  =  p0[e1‘s>1,  eI$2,  ...,  e1$JV]T  and  C  =  cr2Ijv. 

The  elements  of  the  FI  matrix  are 

_  MN  dal  d(Jn  .  2MN  dPo  dPo  .  2  ,fPo  y- 

a £  0Q\  <90„  (j2  dQ\dQu  ^  <90  a  <90„ 


(45) 


(46) 


By  setting  J\ v  =  0,  we  may  determine  the  conditions  underwhich  the  estimates  of  0a  and  0„  will  decouple. 


3.3.1  Plane  Wave 

Let  us  consider  <j>,  9 ,  and  %  as  unknowns,  and  r  and  SNR  as  knowns.  The  estimate  of  \  will  decouple 
from  the  estimates  of  and  9  for  every  value  of  <f>  and  9  if 

N  N  N 

X)  X'i  =  Vi  =  X]  Z'i  =  °’  (47) 

2=1  2=1  2—1 

i.e.,  if  the  center  of  the  array  is  located  at  the  origin.  [If  r  is  unknown,  the  estimate  of  r  will  also  decouple 
from  the  estimates  of  <j>  and  9  if  Eq.  (47)  is  satisfied.]  The  estimates  of  </>  and  9  will  decouple  from  one 
another  if 

N  N  N  N  N 

y'i  =  x'?  and  x'iy'i  =  ^  =  y'^  =  °-  (48) 

2=1  2=1  2=1  2=1  2=1 

Symmetric  planar  array  configurations  such  as  a  circular  array  with  sensors  placed  at  equal  angular  intervals, 
or  a  rectangular  grid  with  sensors  placed  at  the  lattice  points,  meet  these  requirements  provided  that  the 
array  is  located  in  the  ary-plane  and  that  the  array  center  is  located  at  the  origin.  (If  both  r  and  x  are 
unknown,  then  the  quantity  =  +  X  (f^e  phase  of  the  signal  at  the  array  center)  may  be  treated  as 

the  unknown,  and  the  same  results  hold.  Nielsen  [18]  has  performed  an  analysis  for  a  multiple-frequency, 
far-held,  sine  wave  signal  imbedded  in  Gaussian  noise.  The  conditions  he  found  for  the  estimates  of  the 
elevation,  azimuth,  and  phase  at  the  array  center  to  decouple  are  the  same  as  for  the  case  presented  here. 
Among  the  literature  which  examine  array  configurations  that  result  in  the  decoupling  of  the  angle  estimates, 
Refs.  [19,  20,  21]  may  be  of  interest  to  the  reader.) 

Suppose  that  the  center  of  a  sensor  array  is  located  at  the  origin.  Then  if  the  sensors  are  configured  so 
that  Eq.  (48)  is  satisfied, 
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At  9  =  7t/2,  the  CRLB  of  <j>  is  singular.  Therefore,  such  array  configurations  cannot  be  used  to  estimate  the 
azimuthal  AOA  of  a  wave  that  is  propagating  along  the  2:-axis.  Specifically,  a  planar  array  must  be  located 
in  the  xy-plane  in  order  for  Eq.  (48)  to  be  satisfied.  Thus  it  cannot  be  used  to  estimate  the  azimuth  of  a 
normally  incident  wave.  Due  to  the  limitations  of  the  turbulence  model  discussed  in  Sect.  2.1,  we  want  to 
investigate  waves  that  are  nominally  normal  to  a  planar  array.  Therefore  in  the  full  simulation,  the  estimates 
of  the  azimuth  and  elevation  will  always  be  coupled  as  we  cannot  choose  the  xy-plane  as  the  array  plane. 


3.3.2  Spherical  Wave 


Now  let  us  consider  a  spherical  wave.  Suppose  that  r  and  SNR  are  known,  and  that  (j),  0,  and  y  are 
unknown.  The  estimates  of  </>  and  9  will  decouple  from  the  estimate  of  \  f°r  every  value  <f>  and  6  if 
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Therefore,  unlike  the  plane- wave  case,  there  is  no  simple  array  geometry  that  will  result  in  the  decoupling 
of  the  estimates  of  the  AOAs  from  the  estimate  of  the  phase  angle.  The  conditions  for  the  estimates  of  (j) 
and  6  to  decouple  are  also  unattainable  in  practice: 
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It  can  also  be  shown  that  if  r  is  also  unknown,  its  estimate  is  coupled  to  that  of  <j>,  9 ,  and  y. 


3.4  Turbulence 

3.4.1  Plane  Wave 

Let  us  now  consider  a  plane  wave  propagating  through  turbulence.  Unfortunately,  this  case  does  not  lend 
itself  to  being  evaluated  analytically.  We  assume  that  the  wave  is  propagating  near  the  x-axis  and  thus  take 
the  j/2-plane  to  be  the  array  plane.  The  origin  is  taken  to  be  at  the  center  of  the  array.  As  noted  earlier, 
the  estimates  of  9  and  (f>  are  coupled  for  this  choice.  Numerically,  the  estimates  of  9  and  <j>  are  found  to  be 
uncoupled  from  the  estimates  of  l,  c2,  and  SNR.  However  we  find,  again  numerically,  that  placing  the  center 
of  the  array  at  the  origin  is  no  longer  a  sufficient  condition  for  the  estimates  of  <j>  and  9  to  decouple  from  the 
estimates  of  r  and  y:  a  symmetric  array  configuration  such  as  a  rectangular  grid  with  mirror  symmetry  in 
y  about  the  origin  and  mirror  symmetry  in  2  about  the  origin  must  be  used,  though  the  stronger  condition 
of  uniform  spacing  is  not  necessary.  The  results  are  independent  of  the  value  of  y.  (This  is  expected  as  the 
value  of  the  source  phase  should  not  effect  the  estimates  of  the  other  parameters.  Close  inspection  of  the 
second  term  of  Eq.  (40)  reveals  that  its  dependence  in  the  FI  should  cancel.)  The  estimates  of  r,  l,  c2,  and 
SNR  are  all  coupled. 


3.5  Spherical  Wave 

For  the  spherical-wave  analysis,  we  again  assume  that  the  wave  is  propagating  near  the  x-axis,  take  the 
array  plane  to  be  yz-plane,  and  take  the  origin  to  be  at  the  array  center.  Numerically,  we  find:  the  estimates 
of  (j),  9 ,  r  and  y  are  all  coupled;  the  estimates  of  4>  and  9  are  uncoupled  from  the  estimates  of  l,  c2,  and 
SNR0;  and  the  estimates  of  r,  l,  c2,  and  SNR0  are  all  coupled.  Again,  the  results  are  independent  of  the 
value  of  y. 


4  Results 

For  both  the  plane-wave  and  spherical-wave  analyses,  the  array  geometry  considered  is  a  4  x  4  square  grid 
with  spacing  of  d.  In  all  figures,  d/X  =  0.5.  As  the  CRLB  (a)  for  M  independent  and  identically  distributed 
datasets  is  1/yfM  times  the  CRLB  for  one  dataset,  all  results  are  presented  for  M  =  1. 

4.1  Plane  Wave 

The  values  of  cr^,  cr g,  and  ag0  are  the  same  for  normal  incidence  due  to  symmetry.  In  Fig.  6, 

for  normal  incidence  is  plotted  as  a  function  of  c2  and  1/ A  for  r/A  =  500  and  SNR  =  10  dB.  A  vector  von 
Karman  turbulence  spectrum  is  used.  A  peak  is  evident  in  cr^  at  large  c2  and  small  1/ A.  In  this  region 
both  Ti:j  and  rm;n  are  approaching  their  mimimum  values  as  functions  of  the  turbulence  parameters.  In 
fact,  in  the  limit  that  rtJ-  and  rmin  simultaneously  approach  zero,  cr^  — »  oo  as  C  — >  er2I jv  and  /x  — >  0.  For 
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Figure  6:  CRLB  of  < f>  for  a  plane  wave  as  a  function  of  turbulence  parameters.  Calculation  is  for  normal 
incidence,  r/A  =  500,  SNR  =  10  dB,  and  a  vector  von  Karman  spectrum. 


small  values  of  <r2  and  1/ A,  both  F,,  and  rmjn  are  approaching  the  maximum  value  of  1.0,  and  hence  cr^ 
is  approaching  the  limit  for  no  turbulence.  The  behavior  for  other  values  of  the  SNR  is  similar,  with 
increasing  with  decreasing  SNR.  The  corresponding  percent  difference  between  calculated  with  a  zero 
mean  and  with  a  deterministic  mean  is  plotted  in  Fig.  7.  Using  a  deterministic  mean  reduces  the  CRLB. 
The  largest  percent  difference  occurs  for  large  c2  and  small  1/ A.  In  this  region,  however,  the  CRLB  is  large 
and  the  AO  As  cannot  be  estimated  accurately.  In  the  regions  where  rmin  «  1,  the  percent  difference  is 
small,  though  nonzero.  In  the  regions  where  rm;n  ss  0,  the  percent  difference  is  zero,  as  it  should  be. 


Figure  7:  Percent  difference  of  er^  for  a  plane  wave  calculated  with  a  zero  mean  and  with  a  deterministic 
mean.  Calculation  is  for  normal  incidence,  r/A  =  500,  SNR  =  10  dB,  and  a  vector  von  Karman  spectrum. 


Figure  8:  CRLB  of  (j>  for  a  plane  wave  as  a  function  of  normalized  propagation  distance.  All  curves  are  for 
normal  incidence  and  a  vector  von  Karman  spectrum. 


In  Fig.  8,  (T0  for  normal  incidence  is  examined  as  a  function  of  the  normalized  propagation  distance  r/X 
for  a  couple  of  values  of  the  SNR,  <j2,  and  1/ A.  A  vector  von  Karman  spectrum  is  used.  The  points  (lines) 
correspond  to  a  SNR  of  5  dB  (10  dB)  for  every  value  of  r/X.  Though  the  graph  extends  to  smaller  values  of 
r/X  than  are  valid  for  the  turbulence  model,  it  is  useful  to  see  the  limiting  behavior  of  our  model.  [Recall 
the  limits  given  in  Sects.  2. 1-2. 2.  In  particular,  the  approximations  in  Eqs.  (13,  14,  18-20)  must  hold.]  For 
<;2  =  10-6  and  r/X  <  400,  is  limited  by  the  SNR.  As  r/X  increases,  is  driven  primarly  by  the  values 
of  the  turbulence  parameters,  in  particular,  the  index-of-refraction  variance.  The  corresponding  percent 
difference  of  calculated  with  a  zero  mean  and  with  a  deterministic  mean  is  shown  in  Fig.  9. 
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Figure  9:  Percent  difference  of  er^  for  a  plane  wave  calculated  with  a  zero  mean  and  with  a  deterministic 
mean  as  a  function  of  normalized  propagation  distance.  All  curves  are  for  normal  incidence  and  a  vector  von 
Karman  spectrum. 


Figure  10:  CRLB  of  <fi  for  a  spherical  wave  as  a  function  of  turbulence  parameters.  Calculation  is  for  normal 
incidence,  r/A  =  500,  SNRq  =  10  dB  at  TZq/X  =  500,  and  a  vector  von  Karman  spectrum. 


All  plots  presented  have  been  for  normal  incidence.  The  results  for  other  values  of  <j>  and  9  are  similar.  As 
<f>  and  9  increase,  <t^  and  og  increase.  The  angular  dependence  of  and  erg  is  discussed  in  Refs.  [22,  23]  for 
the  same  array  geometry  that  is  considered  here.  Of  particular  interest  is  the  coupling  ^g.  In  Refs.  [22,  23], 
£0 e  is  found  to  be  the  same  as  for  the  no  turbulence  case: 


sin2  <j>  sin2  9 


cos2  9  +  sin2  cf>  sin2  9 


i)292  -  \<j)'l92  +  \(j)29A  . 

o  o 


(52) 


regardless  of  the  the  value  of  r/A,  d/ A,  SNR,  1/ A,  or  c2,  and  regardless  of  the  von  Karman  spectrum  used. 
The  coupling  is  weak:  for  <f>,  9  <  15°,  C,$g  ~  10-3. 

A  vector  von  Karman  spectrum  was  considered  in  Figs.  6-9.  The  analogous  results  using  a  scalar  von 
Karman  spectrum  are  similar.  As  a  function  of  turbulence  parameters,  the  shape  of  the  curves  are  nearly 
identical,  but  the  value  of  <r^,  is  for  the  most  part  smaller.  The  largest  difference  occurs  in  the  regions  where 
the  MCF  for  the  scalar  spectrum  is  appreciably  larger  than  the  MCF  for  a  vector  spectrum  (refer  back  to 
Fig.  3).  More  details  of  the  results  for  a  scalar  spectrum  are  given  in  Refs.  [22,  23].  The  percent  difference  of 
a<f>  using  a  zeromean  and  a  deterministic  mean  is  also  smaller;  in  particular,  the  peak  at  large  <j2  and  small 
1/ A  is  reduced.  Again,  the  use  of  a  deterministic  mean  reduces  the  CRLB. 


4.2  Spherical  Wave 

As  we  are  assuming  spherical  propagation  of  the  wave,  the  signal-to-noise  ratio  is  now  dependent  upon 
r.  Recall  that  the  signal-to-noise  ratio  evaluated  at  r  =  TZ®  is  denoted  SNRo- 

Again,  the  values  of  erg ,  <T0O,  and  ag0  are  the  same  at  for  normal  incidence  due  symmetry.  In  Fig.  10, 
a,/,  for  normal  incidence  is  plotted  versus  1/ A  and  <r2  for  r/A  =  500  and  SNRo  =  10  dB  at  IZo/X  =  500.  A 
vector  von  Karman  spectrum  is  used.  The  overall  values  of  <7^,  are  smaller  than  for  the  plane- wave  case. 
This  is  expected  as  the  values  of  the  MCF  for  a  plane  wave  are  smaller  than  for  a  spherical  (refer  back  to 
Figs.  3-4).  The  corresponding  percent  difference  of  calculated  with  a  zero  mean  and  with  a  deterministic 
mean  is  given  in  Fig.  11.  As  with  the  plane  wave  case,  the  use  of  a  deterministic  mean  reduces  the  CRLBs 
of  the  AOAs.  The  results  are  similar  for  other  values  of  the  SNR. 

In  Fig.  12,  £j0  for  normal  incidence  is  plotted  versus  the  normalized  propagation  distance.  Two  values  of 


Figure  11:  Percent  difference  of  for  a  spherical  wave  calculated  with  zero  mean  and  with  a  deterministic 
mean  as  a  function  of  turbulence  parameters.  Calculation  is  for  normal  incidence,  r/A  =  500,  SNRo  =  10  dB 
at  TZq/X  =  500,  and  a  vector  von  Karman  spectrum. 

<j2,  1/ A,  and  SNR0  evaluated  at  TZq/X  =  500  are  considered.  A  vector  von  Karman  spectrum  is  used.  Again, 
the  graph  is  extended  to  include  smaller  values  of  r/A  than  are  valid  for  the  turbulence  model.  At  small 
values  of  r/A,  we  see  that  is  dependent  upon  the  values  of  the  turbulence  parameters  (particularly  c2) 
and  is  independent  of  the  value  of  SNRo-  Note  the  difference  between  the  outward  spherical  propagation 
and  the  plane  wave  propagation  depicted  in  Fig.  8.  The  corresponding  percent  difference  of  calculated 
with  a  zero  mean  and  with  a  deterministic  mean  is  given  in  Fig.  13. 

Again,  only  the  results  for  normal  incidence  have  been  shown.  For  other  values  of  <f>  and  9 ,  the  results  are 
similiar  and  er^  and  erg  increase  with  increasing  </>  and  9.  At  </>  =  9  =  0,  the  coupling  between  the  estimates 
of  4>  and  9  and  the  estimates  of  r  and  \  is  zero.  At  r/A  =  500  and  l/X  =  10,  we  find  that  for  </>,  9  <  15°: 
C<jjr/\  and  Cer/x  are  on  the  order  of  10-9  for  <j2  =  10-6,  and  on  the  order  of  ICC11  for  <j2  =  10-4;  and 
(gx  are  on  the  order  of  10-9  for  ?2  =  1CK6,  and  on  the  order  of  10-2°  for  ?2  =  1CP4;  and  (^g  is  on  the  order 
of  10~3.  Therefore,  the  coupling  between  the  estimates  of  <f>  and  9  and  estimates  of  x  and  r/A  are  negigible 
for  these  cases.  Moreover,  as  r/A  increases,  these  couplings  rapidly  approach  zero.  Reference  [24]  addresses 
the  angular  dependence  and  couplings  in  more  detail. 

As  with  the  plane  wave  case,  the  results  for  a  scalar  von  Karman  spectrum  are  very  similar  to  those  for 
a  vector  spectrum.  The  CRLBs  of  the  AOAs  are  smaller  for  a  scalar  spectrum.  The  use  of  a  deterministic 
mean  reduces  the  CRLBs  of  the  AOAs,  though  percent  difference  is  not  as  pronounced  as  for  the  vector 
spectrum. 


5  Conclusions 

The  performance  bounds  of  acoustic  arrays  operating  in  atmospheric  turbulence  with  fluctuations  de¬ 
scribed  by  a  von  Karman  spectrum  have  been  examined.  This  analysis  features  four  main  improvements 
upon  earlier  work:  (1)  The  performance  bounds  have  been  generalized  to  weak  as  well  as  strong  scattering. 
(2)  Multiple  unknowns  such  as  the  propagation  distance  of  the  wavefront,  turbulence  parameters,  phase  of 
the  source,  and  signal-to-noise  ratio  have  been  incorporated.  (3)  AOA  estimates  for  three-dimensional  prob¬ 
lems  (i.e.,  two  bearing  angles)  have  been  considered.  (4)  A  multivariate  analysis  for  an  incident  spherical 


Figure  12:  CRLB  of  <j>  for  a  spherical  wave  as  a  function  of  normalized  propagation  distance.  Calculation  is 
for  normal  incidence,  SNRq  evaluated  at  TZq/X  =  500,  and  a  vector  von  Karman  spectrum. 


wave  has  been  developed  in  addition  to  that  for  an  incident  plane  wave. 

The  primary  interest  was  to  analyze  the  Cramer-Rao  lower  bounds  of  the  angles  of  arrival.  For  an  incident 
plane  wave,  we  have  found  that  an  appropriate  choice  of  coordinate  system  and  array  geometry  leads  to 
the  decoupling  of  the  estimates  of  the  AOAs  from  the  estimates  of  the  other  parameters:  the  normalized 
propagation  distance  (r/A),  SNR,  turbulence  parameters  (1/ A  and  <r2),  and  phase  angle  of  the  source  (x).  In 
order  to  remain  consistent  with  the  small-angle  approximation,  we  had  to  choose  a  coordinate  system  that 


Figure  13:  Percent  difference  of  a $  for  a  spherical  wave  calculated  with  a  zero  mean  and  with  a  deterministic 
mean  as  a  function  of  normalized  propagation  distance.  Calculation  is  for  normal  incidence,  SNRo  evaluated 
at  TZq/\  =  500,  and  a  vector  von  Karman  spectrum. 


resulted  in  the  coupling  of  the  estimates  of  the  azimuth  and  zenith;  the  coupling,  however,  was  small.  For 
large  values  of  the  index-of-refraction  variance  and  moderate  to  small  values  of  the  normalized  length  scale, 
we  have  found  that  the  CRLBs  of  the  AOAs  increase  significantly  at  large  propagation  distances.  However, 
for  smaller  values  of  the  index-of-refraction  variance  and  normalized  propagation  distance,  the  SNR  is  the 
limiting  factor. 

For  an  incident  spherical  wave,  we  have  found  that  the  estimates  of  the  AOAs  are  coupled  to  the  estimates 
of  the  normalized  propagation  distance,  the  phase  angle  of  the  source,  and  to  one  another.  The  estimates 
of  the  AOAs  are  uncoupled  from  the  estimates  of  the  turbulence  parameters  and  SNRo-  For  the  array 
geometry  considered,  the  couplings  between  the  estimates  of  the  AOAs  and  the  estimates  of  source  phase 
and  normalized  propagation  distance  are  negigible,  and  the  coupling  of  the  estimates  of  the  azimuth  and 
elevation  is  weak.  For  small  values  of  the  normalized  propagation  distance,  the  CRLBs  of  the  AOAs  are 
dependent  upon  turbulence  parameters  and  independent  of  SNRo- 

We  have  found  that  for  both  a  plane  wave  and  a  spherical  wave,  the  use  of  a  deterministic  mean  reduces 
the  CRLBs  of  the  AOAs  from  that  calculated  with  a  zero  mean.  The  effect  is  more  prevelant  for  a  vector  von 
Karman  spectrum  than  for  a  scalar.  However,  the  largest  percent  difference  occurs  in  the  regions  where  the 
CRLBs  are  large,  and  hence  the  AOAs  cannot  be  well  estimated.  Therefore,  for  this  analysis  the  inclusion  of 
a  deterministic  mean  is  not  very  significant.  However,  for  other  models  of  the  MCF,  the  percent  difference 
may  be  more  significant,  and  thus  a  deterministic  mean  should  be  considered. 

The  results  of  this  analysis  demonstrate  that  scattering  by  atmospheric  turbulence  significantly  affects 
the  performance  of  acoustic  sensor  arrays.  In  order  to  understand  and  circumvent  limitations  on  Army 
acoustical  tracking  systems,  it  is  necessary  to  predict  the  performance  of  acoustic  sensor  arrays  for  various 
atmospheric  conditions.  This  analysis  clearly  demonstrates  the  atmospheric  conditions  that  are  unfavorable 
for  accurate  acoustical  tracking.  While  only  a  single  array  with  a  simple  geometry  was  considered  in  this 
analysis,  the  results  indicate  that  the  effects  of  atmospheric  turbulence  should  be  included  in  performance 
bounds  calculations  for  other  more  complicated  systems  as  well. 

This  analysis  would  benefit  from  an  improved  model  of  the  second  moment  for  oblique  incidence.  Future 
efforts  should  attempt  to  incorporate  the  additional  phenomena  of  ground  reflections  and  refraction  by 
atmospheric  wind  and  temperature  gradients.  These  phenomena  will  likely  have  a  considerable  impact  on 
the  ability  to  estimate  the  elevation.  Numerical  techniques  will  be  required  to  model  these  effects. 
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